Spectral Products¶

Spectral products are created using mathematical combinations of specific bands (wavelength components). These spectral products can be useful for identifying spatial variations in vegetation, water or urbanization. This notebook covers the following spectral products: Fractional Cover, NDBI, NDVI, NDWI, SAVI, EVI

In [1]:
import datacube
dc = datacube.Datacube(app='Spectral_Products')

import sys, os
os.environ['USE_PYGEOS'] = '0'

import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np  
import xarray as xr  

from datacube.utils import masking
from dea_tools.plotting import rgb, display_map
from odc.algo import to_f32

### EASI tools
sys.path.append(os.path.expanduser('../scripts'))
from ceos_utils.data_cube_utilities.clean_mask import landsat_clean_mask_invalid, landsat_qa_clean_mask
from easi_tools import EasiDefaults
from easi_tools import notebook_utils
easi = EasiDefaults() # Get the default parameters for this system
Successfully found configuration for deployment "eail"
In [2]:
cluster, client = notebook_utils.initialize_dask(use_gateway=False)
display(cluster if cluster else client)
print(notebook_utils.localcluster_dashboard(client, server=easi.hub))

LocalCluster

a2da3b55

Dashboard: http://127.0.0.1:8787/status Workers: 4
Total threads: 8 Total memory: 29.00 GiB
Status: running Using processes: True

Scheduler Info

Scheduler

Scheduler-b27d8ebe-785b-4f5d-b2c3-e99981e63dc1

Comm: tcp://127.0.0.1:36521 Workers: 4
Dashboard: http://127.0.0.1:8787/status Total threads: 8
Started: Just now Total memory: 29.00 GiB

Workers

Worker: 0

Comm: tcp://127.0.0.1:36751 Total threads: 2
Dashboard: http://127.0.0.1:39309/status Memory: 7.25 GiB
Nanny: tcp://127.0.0.1:40735
Local directory: /tmp/dask-worker-space/worker-227hjslp

Worker: 1

Comm: tcp://127.0.0.1:35963 Total threads: 2
Dashboard: http://127.0.0.1:42569/status Memory: 7.25 GiB
Nanny: tcp://127.0.0.1:40263
Local directory: /tmp/dask-worker-space/worker-9raofrjh

Worker: 2

Comm: tcp://127.0.0.1:36497 Total threads: 2
Dashboard: http://127.0.0.1:42951/status Memory: 7.25 GiB
Nanny: tcp://127.0.0.1:46399
Local directory: /tmp/dask-worker-space/worker-pfff9htp

Worker: 3

Comm: tcp://127.0.0.1:35369 Total threads: 2
Dashboard: http://127.0.0.1:34953/status Memory: 7.25 GiB
Nanny: tcp://127.0.0.1:36405
Local directory: /tmp/dask-worker-space/worker-hrwxic0a
https://hub.eail.easi-eo.solutions/user/jhodge/proxy/8787/status
In [3]:
from datacube.utils.aws import configure_s3_access
configure_s3_access(aws_unsigned=False, requester_pays=True, client=client)
Out[3]:
<botocore.credentials.Credentials at 0x7fd232798340>
In [4]:
# Select a Product
product = "landsat8_c2l2_sr"
In [5]:
# Select an analysis region (Lat-Lon) within the extents listed above. 
# Select a time period (Min-Max) within the extents listed above (Year-Month-Day)
# This region and time period will be used for the cloud assessment

# Kigoma, Tanzania
# latitude = (-4.95, -4.84) 
# longitude = (29.59, 29.75) 

# Dar es Salaam, Tanzania
# latitude = (-7.0, -6.75)
# longitude = (39.1, 39.35)

# Mtera Reservoir, Tanzania
# latitude = (-7.22, -6.80) 
# longitude = (35.60, 36.00) 

# Mining Region near Obuasi, Ghana
# latitude = (6.0985, 6.2675)
# longitude = (-2.050, -1.8629)

# Efate Island, Vanuatu (near Port Vila)
# latitude = (-17.810, -17.691) 
# longitude = (168.247, 168.349) 

# Port Vila, Vanuatu
# latitude = (-17.8, -17.65) 
# longitude = (168.25, 168.375) 

# Solomons Island - Coconut Crop Damaged by Beetles - FULL
# latitude = (-9.4412, -9.4357)
# longitude = (160.0241, 160.0288)

# Solomons Island - Coconut Crop Damaged by Beetles - MIDDLE
# latitude = (-9.4391, -9.4372)
# longitude = (160.0257, 160.0278)

latitude = easi.latitude
longitude = easi.longitude

# Time Period (Year-Mon-Day)
time_extents = ('2015-1-1', '2015-12-31')
In [6]:
# The code below renders a map that can be used to view the region.
display_map(longitude,latitude)
/env/lib/python3.10/site-packages/dea_tools/plotting.py:313: FutureWarning: This function is deprecated. See: https://pyproj4.github.io/pyproj/stable/gotchas.html#upgrading-to-pyproj-2-from-pyproj-1
  all_longitude, all_latitude = transform(Proj(crs), Proj('EPSG:4326'), all_x,
Out[6]:
Make this Notebook Trusted to load map: File -> Trust Notebook

Load the dataset and the required spectral bands or other parameters¶

In [7]:
landsat_dataset = dc.load(latitude = latitude,
                          longitude = longitude,
                          time = time_extents,
                          product = product,
                          measurements = ['red', 'green', 'blue', 'nir', 'swir1', 'swir2', 'pixel_qa'],
                          output_crs = 'EPSG:6933',
                          resolution = (-30,30),
                          dask_chunks = {'time':1}
                         ) 
In [8]:
# Displays the first few values of each data array to check the content
# Latitude and Longitude numbers = number of pixels in each dimension
# Time = number of time slices in the dataset

landsat_dataset
Out[8]:
<xarray.Dataset>
Dimensions:      (time: 36, y: 342, x: 322)
Coordinates:
  * time         (time) datetime64[ns] 2015-01-05T15:46:53.790263 ... 2015-12...
  * y            (y) float64 4.418e+06 4.418e+06 ... 4.408e+06 4.408e+06
  * x            (x) float64 -7.386e+06 -7.386e+06 ... -7.376e+06 -7.376e+06
    spatial_ref  int32 6933
Data variables:
    red          (time, y, x) uint16 dask.array<chunksize=(1, 342, 322), meta=np.ndarray>
    green        (time, y, x) uint16 dask.array<chunksize=(1, 342, 322), meta=np.ndarray>
    blue         (time, y, x) uint16 dask.array<chunksize=(1, 342, 322), meta=np.ndarray>
    nir          (time, y, x) uint16 dask.array<chunksize=(1, 342, 322), meta=np.ndarray>
    swir1        (time, y, x) uint16 dask.array<chunksize=(1, 342, 322), meta=np.ndarray>
    swir2        (time, y, x) uint16 dask.array<chunksize=(1, 342, 322), meta=np.ndarray>
    pixel_qa     (time, y, x) uint16 dask.array<chunksize=(1, 342, 322), meta=np.ndarray>
Attributes:
    crs:           EPSG:6933
    grid_mapping:  spatial_ref
xarray.Dataset
    • time: 36
    • y: 342
    • x: 322
    • time
      (time)
      datetime64[ns]
      2015-01-05T15:46:53.790263 ... 2...
      units :
      seconds since 1970-01-01 00:00:00
      array(['2015-01-05T15:46:53.790263000', '2015-01-30T15:40:37.355183000',
             '2015-02-06T15:46:45.651144000', '2015-02-15T15:40:28.526601000',
             '2015-03-19T15:40:14.706886000', '2015-03-26T15:46:21.170165000',
             '2015-04-04T15:40:02.713256000', '2015-04-11T15:46:13.625587000',
             '2015-04-20T15:40:02.235599000', '2015-04-27T15:46:08.223358000',
             '2015-05-06T15:39:47.603585000', '2015-05-13T15:45:52.093905000',
             '2015-05-22T15:39:43.391309000', '2015-05-29T15:45:55.414220000',
             '2015-06-07T15:39:52.787430000', '2015-06-14T15:46:07.814585000',
             '2015-06-23T15:39:59.338668000', '2015-06-30T15:46:14.742711000',
             '2015-07-09T15:40:09.955667000', '2015-07-16T15:46:24.314629000',
             '2015-07-25T15:40:16.195856000', '2015-08-01T15:46:28.107825000',
             '2015-08-17T15:46:34.758025000', '2015-08-26T15:40:27.437353000',
             '2015-09-02T15:46:39.676433000', '2015-09-11T15:40:34.119112000',
             '2015-09-18T15:46:48.349171000', '2015-09-27T15:40:40.387129000',
             '2015-10-13T15:40:40.514877000', '2015-10-20T15:46:54.010262000',
             '2015-10-29T15:40:46.414136000', '2015-11-05T15:46:57.863683000',
             '2015-11-14T15:40:46.710019000', '2015-11-21T15:46:59.672526000',
             '2015-12-07T15:46:58.456582000', '2015-12-16T15:40:48.214294000'],
            dtype='datetime64[ns]')
    • y
      (y)
      float64
      4.418e+06 4.418e+06 ... 4.408e+06
      units :
      metre
      resolution :
      -30.0
      crs :
      EPSG:6933
      array([4418325., 4418295., 4418265., ..., 4408155., 4408125., 4408095.])
    • x
      (x)
      float64
      -7.386e+06 ... -7.376e+06
      units :
      metre
      resolution :
      30.0
      crs :
      EPSG:6933
      array([-7386015., -7385985., -7385955., ..., -7376445., -7376415., -7376385.])
    • spatial_ref
      ()
      int32
      6933
      spatial_ref :
      PROJCS["WGS 84 / NSIDC EASE-Grid 2.0 Global",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Cylindrical_Equal_Area"],PARAMETER["standard_parallel_1",30],PARAMETER["central_meridian",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","6933"]]
      grid_mapping_name :
      lambert_cylindrical_equal_area
      array(6933, dtype=int32)
    • red
      (time, y, x)
      uint16
      dask.array<chunksize=(1, 342, 322), meta=np.ndarray>
      units :
      reflectance
      nodata :
      0
      crs :
      EPSG:6933
      grid_mapping :
      spatial_ref
      Array Chunk
      Bytes 7.56 MiB 215.09 kiB
      Shape (36, 342, 322) (1, 342, 322)
      Dask graph 36 chunks in 1 graph layer
      Data type uint16 numpy.ndarray
      322 342 36
    • green
      (time, y, x)
      uint16
      dask.array<chunksize=(1, 342, 322), meta=np.ndarray>
      units :
      reflectance
      nodata :
      0
      crs :
      EPSG:6933
      grid_mapping :
      spatial_ref
      Array Chunk
      Bytes 7.56 MiB 215.09 kiB
      Shape (36, 342, 322) (1, 342, 322)
      Dask graph 36 chunks in 1 graph layer
      Data type uint16 numpy.ndarray
      322 342 36
    • blue
      (time, y, x)
      uint16
      dask.array<chunksize=(1, 342, 322), meta=np.ndarray>
      units :
      reflectance
      nodata :
      0
      crs :
      EPSG:6933
      grid_mapping :
      spatial_ref
      Array Chunk
      Bytes 7.56 MiB 215.09 kiB
      Shape (36, 342, 322) (1, 342, 322)
      Dask graph 36 chunks in 1 graph layer
      Data type uint16 numpy.ndarray
      322 342 36
    • nir
      (time, y, x)
      uint16
      dask.array<chunksize=(1, 342, 322), meta=np.ndarray>
      units :
      reflectance
      nodata :
      0
      crs :
      EPSG:6933
      grid_mapping :
      spatial_ref
      Array Chunk
      Bytes 7.56 MiB 215.09 kiB
      Shape (36, 342, 322) (1, 342, 322)
      Dask graph 36 chunks in 1 graph layer
      Data type uint16 numpy.ndarray
      322 342 36
    • swir1
      (time, y, x)
      uint16
      dask.array<chunksize=(1, 342, 322), meta=np.ndarray>
      units :
      reflectance
      nodata :
      0
      crs :
      EPSG:6933
      grid_mapping :
      spatial_ref
      Array Chunk
      Bytes 7.56 MiB 215.09 kiB
      Shape (36, 342, 322) (1, 342, 322)
      Dask graph 36 chunks in 1 graph layer
      Data type uint16 numpy.ndarray
      322 342 36
    • swir2
      (time, y, x)
      uint16
      dask.array<chunksize=(1, 342, 322), meta=np.ndarray>
      units :
      reflectance
      nodata :
      0
      crs :
      EPSG:6933
      grid_mapping :
      spatial_ref
      Array Chunk
      Bytes 7.56 MiB 215.09 kiB
      Shape (36, 342, 322) (1, 342, 322)
      Dask graph 36 chunks in 1 graph layer
      Data type uint16 numpy.ndarray
      322 342 36
    • pixel_qa
      (time, y, x)
      uint16
      dask.array<chunksize=(1, 342, 322), meta=np.ndarray>
      units :
      bit_index
      nodata :
      1
      flags_definition :
      {'snow': {'bits': 5, 'values': {'0': 'not_high_confidence', '1': 'high_confidence'}}, 'clear': {'bits': 6, 'values': {'0': 'not_clear', '1': 'clear'}}, 'cloud': {'bits': 3, 'values': {'0': 'not_high_confidence', '1': 'high_confidence'}}, 'water': {'bits': 7, 'values': {'0': 'land_or_cloud', '1': 'water'}}, 'cirrus': {'bits': 2, 'values': {'0': 'not_high_confidence', '1': 'high_confidence'}}, 'nodata': {'bits': 0, 'values': {'0': False, '1': True}}, 'qa_pixel': {'bits': [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15], 'values': {'1': 'Fill', '2': 'Dilated Cloud', '4': 'Cirrus', '8': 'Cloud', '16': 'Cloud Shadow', '32': 'Snow', '64': 'Clear', '128': 'Water', '256': 'Cloud Confidence low bit', '512': 'Cloud Confidence high bit', '1024': 'Cloud Shadow Confidence low bit', '2048': 'Cloud Shadow Confidence high bit', '4096': 'Snow Ice Confidence low bit', '8192': 'Snow Ice Confidence high bit', '16384': 'Cirrus Confidence low bit', '32768': 'Cirrus Confidence high bit'}, 'description': 'Level 2 pixel quality'}, 'cloud_shadow': {'bits': 4, 'values': {'0': 'not_high_confidence', '1': 'high_confidence'}}, 'dilated_cloud': {'bits': 1, 'values': {'0': 'not_dilated', '1': 'dilated'}}, 'cloud_confidence': {'bits': [8, 9], 'values': {'0': 'none', '1': 'low', '2': 'medium', '3': 'high'}}, 'cirrus_confidence': {'bits': [14, 15], 'values': {'0': 'none', '1': 'low', '2': 'reserved', '3': 'high'}}, 'snow_ice_confidence': {'bits': [12, 13], 'values': {'0': 'none', '1': 'low', '2': 'reserved', '3': 'high'}}, 'cloud_shadow_confidence': {'bits': [10, 11], 'values': {'0': 'none', '1': 'low', '2': 'reserved', '3': 'high'}}}
      crs :
      EPSG:6933
      grid_mapping :
      spatial_ref
      Array Chunk
      Bytes 7.56 MiB 215.09 kiB
      Shape (36, 342, 322) (1, 342, 322)
      Dask graph 36 chunks in 1 graph layer
      Data type uint16 numpy.ndarray
      322 342 36
    • time
      PandasIndex
      PandasIndex(DatetimeIndex(['2015-01-05 15:46:53.790263', '2015-01-30 15:40:37.355183',
                     '2015-02-06 15:46:45.651144', '2015-02-15 15:40:28.526601',
                     '2015-03-19 15:40:14.706886', '2015-03-26 15:46:21.170165',
                     '2015-04-04 15:40:02.713256', '2015-04-11 15:46:13.625587',
                     '2015-04-20 15:40:02.235599', '2015-04-27 15:46:08.223358',
                     '2015-05-06 15:39:47.603585', '2015-05-13 15:45:52.093905',
                     '2015-05-22 15:39:43.391309', '2015-05-29 15:45:55.414220',
                     '2015-06-07 15:39:52.787430', '2015-06-14 15:46:07.814585',
                     '2015-06-23 15:39:59.338668', '2015-06-30 15:46:14.742711',
                     '2015-07-09 15:40:09.955667', '2015-07-16 15:46:24.314629',
                     '2015-07-25 15:40:16.195856', '2015-08-01 15:46:28.107825',
                     '2015-08-17 15:46:34.758025', '2015-08-26 15:40:27.437353',
                     '2015-09-02 15:46:39.676433', '2015-09-11 15:40:34.119112',
                     '2015-09-18 15:46:48.349171', '2015-09-27 15:40:40.387129',
                     '2015-10-13 15:40:40.514877', '2015-10-20 15:46:54.010262',
                     '2015-10-29 15:40:46.414136', '2015-11-05 15:46:57.863683',
                     '2015-11-14 15:40:46.710019', '2015-11-21 15:46:59.672526',
                     '2015-12-07 15:46:58.456582', '2015-12-16 15:40:48.214294'],
                    dtype='datetime64[ns]', name='time', freq=None))
    • y
      PandasIndex
      PandasIndex(Float64Index([4418325.0, 4418295.0, 4418265.0, 4418235.0, 4418205.0, 4418175.0,
                    4418145.0, 4418115.0, 4418085.0, 4418055.0,
                    ...
                    4408365.0, 4408335.0, 4408305.0, 4408275.0, 4408245.0, 4408215.0,
                    4408185.0, 4408155.0, 4408125.0, 4408095.0],
                   dtype='float64', name='y', length=342))
    • x
      PandasIndex
      PandasIndex(Float64Index([-7386015.0, -7385985.0, -7385955.0, -7385925.0, -7385895.0,
                    -7385865.0, -7385835.0, -7385805.0, -7385775.0, -7385745.0,
                    ...
                    -7376655.0, -7376625.0, -7376595.0, -7376565.0, -7376535.0,
                    -7376505.0, -7376475.0, -7376445.0, -7376415.0, -7376385.0],
                   dtype='float64', name='x', length=322))
  • crs :
    EPSG:6933
    grid_mapping :
    spatial_ref

Mask out clouds and cloud shadows + water (if desired) and create a median mosaic¶

In [9]:
cloud_mask = masking.make_mask(landsat_dataset['pixel_qa'], cloud='not_high_confidence', cloud_shadow='not_high_confidence')

land_mask = masking.make_mask(landsat_dataset['pixel_qa'], water='land_or_cloud')

# Land and Water Dataset = Land and Water pixels with NO Clouds and NO Cloud Shadows
land_and_water_dataset = landsat_dataset.where(cloud_mask)

# Land Dataset = Land ONLY pixels with NO Clouds, NO Cloud Shadows and NO Water pixels
land_dataset = landsat_dataset.where(land_mask)

from ceos_utils.data_cube_utilities.dc_mosaic import create_median_mosaic, create_max_ndvi_mosaic, create_hdmedians_multiple_band_mosaic
In [10]:
# Select a compositing method to create your cloud-filtered mosaic
# Remove the comments from the pair of lines under one of the mosaic types
# Options are: Median or Max_NDVI 

# This is the MEDIAN mosaic
scale = 0.0275
offset = 0.2
land_and_water_composite = to_f32(create_median_mosaic(land_and_water_dataset, cloud_mask), scale=scale, offset=offset)
land_composite = to_f32(create_median_mosaic(land_dataset, land_mask), scale=scale, offset=offset)
cloud_mask_composite = cloud_mask.max('time')

# This is the MAX_NDVI mosaic
# land_and_water_composite = create_max_ndvi_mosaic(land_and_water_dataset, cloud_mask)
# land_composite = create_max_ndvi_mosaic(land_dataset, land_mask)
In [11]:
# Show the land-only composite (water will be removed ... black pixels)

# RGB image options
# Standard RGB = 321 = Red, Green, Blue
# False Color = 543 = SWIR1, NIR, Red
# False Color (Landsat Mosaic) = 742 = SWIR2, NIR, Green

rgb(land_composite, bands=['swir2', 'nir', 'green'], size=8)
plt.axis('off')
plt.show()

Land Spectral Indices¶

In [12]:
def NDBI(dataset):
        return (dataset.swir1 - dataset.nir)/(dataset.swir1 + dataset.nir)
In [13]:
def NDVI(dataset):
    return (dataset.nir - dataset.red)/(dataset.nir + dataset.red)
In [14]:
def NDWI(dataset):
    return (dataset.green - dataset.nir)/(dataset.green + dataset.nir)
In [15]:
def SAVI(dataset):
        return (dataset.nir - dataset.red)/(dataset.nir + dataset.red + 0.5)*1.5
In [16]:
def EVI(dataset):
        return 2.5*(dataset.nir - dataset.red)/(dataset.nir + 6.0*dataset.red - 7.5*dataset.blue + 1.0)
In [17]:
# Water pixels masked (land only)
ndbi = NDBI(land_composite)  # Normalized Difference Build Up (Urbanization) Index
ndvi = NDVI(land_composite)  # Normalized Difference Vegetation Index
ndwi = NDWI(land_composite) # Normalized Difference Water Index

# Water pixels not masked
ndbi2 = NDBI(land_and_water_composite)  # Normalized Difference Build Up (Urbanization) Index
ndvi2 = NDVI(land_and_water_composite)  # Normalized Difference Vegetation Index
ndwi2 = NDWI(land_and_water_composite) # Normalized Difference Water Index

# Other Indices
savi = SAVI(land_composite)  # Soil Adjusted Vegetation Index 
evi = EVI(land_composite) # Enhanced Vegetation Index
In [18]:
ds_ndvi = ndvi2.to_dataset(name = "NDVI")
ds_ndwi = ndwi2.to_dataset(name=  "NDWI")
ds_ndbi = ndbi2.to_dataset(name = "NDBI")
normalization_dataset = ds_ndvi.merge(ds_ndwi).merge(ds_ndbi)

Land Fractional Cover¶

Fractional Cover (FC) is used for landcover type estimation (vegetation, non-green vegetation, bare soil) of each pixel. We use a model from CSIRO (Juan Gerschmann) and apply it to a median mosaic where: Bare Soil = bs, Photosynthetic Vegetation = pv, and Non-Photosynthetic Vegetation = npv. The product is a False Color RGB result where RGB = bs/pv/npv.

In [19]:
land_composite = land_composite.rename_dims({'x':'longitude','y':'latitude'})
In [20]:
land_composite
Out[20]:
<xarray.Dataset>
Dimensions:      (latitude: 342, longitude: 322)
Coordinates:
  * y            (latitude) float64 4.418e+06 4.418e+06 ... 4.408e+06 4.408e+06
  * x            (longitude) float64 -7.386e+06 -7.386e+06 ... -7.376e+06
    spatial_ref  int32 6933
Dimensions without coordinates: latitude, longitude
Data variables:
    red          (latitude, longitude) float32 dask.array<chunksize=(342, 322), meta=np.ndarray>
    green        (latitude, longitude) float32 dask.array<chunksize=(342, 322), meta=np.ndarray>
    blue         (latitude, longitude) float32 dask.array<chunksize=(342, 322), meta=np.ndarray>
    nir          (latitude, longitude) float32 dask.array<chunksize=(342, 322), meta=np.ndarray>
    swir1        (latitude, longitude) float32 dask.array<chunksize=(342, 322), meta=np.ndarray>
    swir2        (latitude, longitude) float32 dask.array<chunksize=(342, 322), meta=np.ndarray>
    pixel_qa     (latitude, longitude) float32 dask.array<chunksize=(342, 322), meta=np.ndarray>
xarray.Dataset
    • latitude: 342
    • longitude: 322
    • y
      (latitude)
      float64
      4.418e+06 4.418e+06 ... 4.408e+06
      units :
      metre
      resolution :
      -30.0
      crs :
      EPSG:6933
      array([4418325., 4418295., 4418265., ..., 4408155., 4408125., 4408095.])
    • x
      (longitude)
      float64
      -7.386e+06 ... -7.376e+06
      units :
      metre
      resolution :
      30.0
      crs :
      EPSG:6933
      array([-7386015., -7385985., -7385955., ..., -7376445., -7376415., -7376385.])
    • spatial_ref
      ()
      int32
      6933
      spatial_ref :
      PROJCS["WGS 84 / NSIDC EASE-Grid 2.0 Global",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Cylindrical_Equal_Area"],PARAMETER["standard_parallel_1",30],PARAMETER["central_meridian",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","6933"]]
      grid_mapping_name :
      lambert_cylindrical_equal_area
      array(6933, dtype=int32)
    • red
      (latitude, longitude)
      float32
      dask.array<chunksize=(342, 322), meta=np.ndarray>
      Array Chunk
      Bytes 430.17 kiB 430.17 kiB
      Shape (342, 322) (342, 322)
      Dask graph 1 chunks in 17 graph layers
      Data type float32 numpy.ndarray
      322 342
    • green
      (latitude, longitude)
      float32
      dask.array<chunksize=(342, 322), meta=np.ndarray>
      Array Chunk
      Bytes 430.17 kiB 430.17 kiB
      Shape (342, 322) (342, 322)
      Dask graph 1 chunks in 17 graph layers
      Data type float32 numpy.ndarray
      322 342
    • blue
      (latitude, longitude)
      float32
      dask.array<chunksize=(342, 322), meta=np.ndarray>
      Array Chunk
      Bytes 430.17 kiB 430.17 kiB
      Shape (342, 322) (342, 322)
      Dask graph 1 chunks in 17 graph layers
      Data type float32 numpy.ndarray
      322 342
    • nir
      (latitude, longitude)
      float32
      dask.array<chunksize=(342, 322), meta=np.ndarray>
      Array Chunk
      Bytes 430.17 kiB 430.17 kiB
      Shape (342, 322) (342, 322)
      Dask graph 1 chunks in 17 graph layers
      Data type float32 numpy.ndarray
      322 342
    • swir1
      (latitude, longitude)
      float32
      dask.array<chunksize=(342, 322), meta=np.ndarray>
      Array Chunk
      Bytes 430.17 kiB 430.17 kiB
      Shape (342, 322) (342, 322)
      Dask graph 1 chunks in 17 graph layers
      Data type float32 numpy.ndarray
      322 342
    • swir2
      (latitude, longitude)
      float32
      dask.array<chunksize=(342, 322), meta=np.ndarray>
      Array Chunk
      Bytes 430.17 kiB 430.17 kiB
      Shape (342, 322) (342, 322)
      Dask graph 1 chunks in 17 graph layers
      Data type float32 numpy.ndarray
      322 342
    • pixel_qa
      (latitude, longitude)
      float32
      dask.array<chunksize=(342, 322), meta=np.ndarray>
      Array Chunk
      Bytes 430.17 kiB 430.17 kiB
      Shape (342, 322) (342, 322)
      Dask graph 1 chunks in 16 graph layers
      Data type float32 numpy.ndarray
      322 342
    • y
      PandasIndex
      PandasIndex(Float64Index([4418325.0, 4418295.0, 4418265.0, 4418235.0, 4418205.0, 4418175.0,
                    4418145.0, 4418115.0, 4418085.0, 4418055.0,
                    ...
                    4408365.0, 4408335.0, 4408305.0, 4408275.0, 4408245.0, 4408215.0,
                    4408185.0, 4408155.0, 4408125.0, 4408095.0],
                   dtype='float64', name='y', length=342))
    • x
      PandasIndex
      PandasIndex(Float64Index([-7386015.0, -7385985.0, -7385955.0, -7385925.0, -7385895.0,
                    -7385865.0, -7385835.0, -7385805.0, -7385775.0, -7385745.0,
                    ...
                    -7376655.0, -7376625.0, -7376595.0, -7376565.0, -7376535.0,
                    -7376505.0, -7376475.0, -7376445.0, -7376415.0, -7376385.0],
                   dtype='float64', name='x', length=322))
NOTE: for some reason this doesn't work properly and the Fractional cover image gets messed up (upside down). There is something going wrong with the handling of arrays. I think it would be better if the function worked in pure xarray rather than numpy arrays as there are too many possiblities for losing the relationships with dimensions.
In [21]:
from ceos_utils.data_cube_utilities.dc_fractional_coverage_classifier import frac_coverage_classify 
frac_classes = frac_coverage_classify(land_composite, clean_mask = 
                                      np.ones(land_composite.pixel_qa.shape).astype(np.bool))
/tmp/ipykernel_472/2442966445.py:3: DeprecationWarning: `np.bool` is a deprecated alias for the builtin `bool`. To silence this warning, use `bool` by itself. Doing this will not modify any behavior and is safe. If you specifically wanted the numpy scalar type, use `np.bool_` here.
Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations
  np.ones(land_composite.pixel_qa.shape).astype(np.bool))
/env/lib/python3.10/site-packages/rasterio/warp.py:344: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned.
  _reproject(
In [22]:
frac_classes
Out[22]:
<xarray.Dataset>
Dimensions:    (latitude: 342, longitude: 322)
Coordinates:
  * latitude   (latitude) int64 0 1 2 3 4 5 6 7 ... 335 336 337 338 339 340 341
  * longitude  (longitude) int64 0 1 2 3 4 5 6 7 ... 315 316 317 318 319 320 321
Data variables:
    bs         (latitude, longitude) float32 67.0 60.0 55.0 ... 62.0 63.0 77.0
    pv         (latitude, longitude) float32 22.0 28.0 32.0 ... 29.0 27.0 14.0
    npv        (latitude, longitude) float32 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
xarray.Dataset
    • latitude: 342
    • longitude: 322
    • latitude
      (latitude)
      int64
      0 1 2 3 4 5 ... 337 338 339 340 341
      array([  0,   1,   2, ..., 339, 340, 341])
    • longitude
      (longitude)
      int64
      0 1 2 3 4 5 ... 317 318 319 320 321
      array([  0,   1,   2, ..., 319, 320, 321])
    • bs
      (latitude, longitude)
      float32
      67.0 60.0 55.0 ... 62.0 63.0 77.0
      array([[67., 60., 55., ..., 73., 66., 62.],
             [62., 66., 64., ..., 68., 65., 69.],
             [63., 65., 64., ..., 79., 79., 79.],
             ...,
             [93., 94., 97., ..., 81., 82., 93.],
             [96., 96., 96., ..., 77., 67., 88.],
             [96., 97., 97., ..., 62., 63., 77.]], dtype=float32)
    • pv
      (latitude, longitude)
      float32
      22.0 28.0 32.0 ... 29.0 27.0 14.0
      array([[22., 28., 32., ..., 17., 24., 27.],
             [26., 22., 24., ..., 23., 25., 21.],
             [26., 24., 25., ..., 12., 12., 12.],
             ...,
             [ 3.,  3.,  0., ..., 12., 10.,  0.],
             [ 0.,  0.,  0., ..., 14., 22.,  4.],
             [ 0.,  0.,  0., ..., 29., 27., 14.]], dtype=float32)
    • npv
      (latitude, longitude)
      float32
      0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
      array([[0., 0., 0., ..., 0., 0., 0.],
             [0., 0., 0., ..., 0., 0., 0.],
             [0., 0., 0., ..., 0., 0., 0.],
             ...,
             [0., 0., 0., ..., 0., 0., 0.],
             [0., 0., 0., ..., 0., 0., 0.],
             [0., 0., 0., ..., 0., 0., 0.]], dtype=float32)
    • latitude
      PandasIndex
      PandasIndex(RangeIndex(start=0, stop=342, step=1, name='latitude'))
    • longitude
      PandasIndex
      PandasIndex(RangeIndex(start=0, stop=322, step=1, name='longitude'))
In [23]:
# Plot of Fractional Cover
# RED = Bare Soil or Urban Areas
# BLUE = Non-Green Vegetation
# GREEN = Green Vegetation
# BLACK = Water

# Plot of RGB = NDBI-NDVI-NDWI
# RED = Bare Soil or Urban Areas
# GREEN = Vegetation
# BLUE = Water

fig, ax = plt.subplots(1, 2, figsize=(16, 8))

fc_rgb = frac_classes[['bs', 'pv', 'npv']].to_array()
ndi_rgb = normalization_dataset[['NDBI','NDVI','NDWI']].to_array()

(fc_rgb).plot.imshow(ax=ax[0], vmin=0.0, vmax=100.0)
(ndi_rgb).plot.imshow(ax=ax[1], vmin=-1.0, vmax=1.0)

ax[0].set_title('Fractional Cover'), ax[0].xaxis.set_visible(False), ax[0].yaxis.set_visible(False)
ax[1].set_title('Normalized Difference RGB'), ax[1].xaxis.set_visible(False), ax[1].yaxis.set_visible(False)
plt.show()
/env/lib/python3.10/site-packages/rasterio/warp.py:344: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned.
  _reproject(
In [24]:
import matplotlib as mpl
# %pylab inline
# matplotlib.rcParams['figure.figsize'] = (10, 10)
In [25]:
# Create a custom colour map for NDVI
# Water (blue) = NDVI -1.0 to 0.05
# Urban or Bare Soil (brown) = NDVI 0.05 to 0.25
# Low Vegetation (tan) = NDVI 0.25 to 0.4
# Croplands (light green) = NDVI 0.4 to 0.6
# Dense Vegetation / Forests (dark green) = NDVI 0.6 to 1.0

ndvi_cmap = mpl.colors.ListedColormap(['blue', '#a52a2a','#ffffcc' ,  '#2eb82e',  '#006600'])
ndvi_bounds = [-1, 0.05, 0.25,  0.4,  0.6, 1]
ndvi_norm = mpl.colors.BoundaryNorm(ndvi_bounds, ndvi_cmap.N)
In [26]:
fig, ax = plt.subplots(1, 2, figsize=(16, 8))
ax[0].imshow(ndvi2, cmap="Greens", vmin=0.0, vmax=1.0)
ax[1].imshow(ndvi2, cmap=ndvi_cmap, norm = ndvi_norm)
ax[0].set_title('NDVI Basic'), ax[0].xaxis.set_visible(False), ax[0].yaxis.set_visible(False)
ax[1].set_title('NDVI Custom Colormap'), ax[1].xaxis.set_visible(False), ax[1].yaxis.set_visible(False)
plt.show()
/env/lib/python3.10/site-packages/rasterio/warp.py:344: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned.
  _reproject(
In [27]:
# EVI and SAVI indices using "land only" pixels
fig, ax = plt.subplots(1, 2, figsize=(18, 8))
(evi).plot.imshow(ax=ax[0], cmap="Greens", vmin=-1.0, vmax=2.5)
(savi).plot.imshow(ax=ax[1], cmap="Greens", vmin=-1.0, vmax=1.0)
ax[0].set_title('EVI'), ax[0].xaxis.set_visible(False), ax[0].yaxis.set_visible(False)
ax[1].set_title('SAVI'), ax[1].xaxis.set_visible(False), ax[1].yaxis.set_visible(False)
plt.show()

Create a threshold plot¶

First we will define a minimum threshold and a maximum threshold. Then you will create a plot that colors the region between the threshold a single color (e.g. red) and the region outside the threshold will be BLACK or WHITE. Also, we will calculate the % of pixels and the number of pixels in the threshold range.

In [28]:
from matplotlib.ticker import FuncFormatter

def threshold_plot(da, min_threshold, max_threshold, mask = None, width = 10, *args, **kwargs): 
    color_in    = np.array([255,0,0])
    color_out   = np.array([0,0,0])
    color_cloud = np.array([255,255,255])
    
    array = np.zeros((*da.values.shape, 3)).astype(np.int16)
    
    inside  = np.logical_and(da.values > min_threshold, da.values < max_threshold)
    outside = np.invert(inside)
    masked  = np.zeros(da.values.shape).astype(bool) if mask is None else mask
    
    array[inside] =  color_in
    array[outside] = color_out
    array[masked] =  color_cloud

    def figure_ratio(ds, fixed_width = 10):
        width = fixed_width
        height = len(ds.latitude) * (fixed_width / len(ds.longitude))
        return (width, height)


    fig, ax = plt.subplots(figsize = figure_ratio(da,fixed_width = width))
    
    lat_formatter = FuncFormatter(lambda y_val, tick_pos: "{0:.3f}".format(da.latitude.values[tick_pos] ))
    lon_formatter = FuncFormatter(lambda x_val, tick_pos: "{0:.3f}".format(da.longitude.values[tick_pos]))

    ax.xaxis.set_major_formatter(lon_formatter)
    ax.yaxis.set_major_formatter(lat_formatter)
    
    plt.title("Threshold: {} < x < {}".format(min_threshold, max_threshold))
    plt.xlabel('Longitude')
    plt.ylabel('Latitude')
    
    plt.imshow(array, *args, **kwargs)
    plt.axis('off')
    plt.show()
In [29]:
# Select a threshold range for your spectral variable and generate a plot
# Remove comments from the set of 3 lines for your desired variable
# Variable choices are NDBI, NDVI, EVI, FC-Bare Soil, FC-Photosynthetic Vegetation

# NDBI (Buildup Index) = -1.0 to 1.0 (full range)
# NDBI 0.0 to 0.2 is typical for urban areas
# -----------------------
# minimum_threshold = 0.0
# maximum_threshold = 0.3
# threshold_plot(ndbi, minimum_threshold, maximum_threshold, width = 8)

# NDVI (Vegetation Index) = -1.0 to 1.0
# NDVI < 0.0 = non-vegetation (bare soil)
# NDVI 0.2 to 0.6 = grasslands
# NDVI 0.6 to 0.9 = dense vegetation / trees
# -----------------------
minimum_threshold = 0.05
maximum_threshold = 0.25
threshold_plot(ndvi.rename({'x':'longitude','y':'latitude'}), minimum_threshold, maximum_threshold, width = 8)

# EVI (Vegetation Index) = -1.0 to 2.5
# EVI 2.0 to 2.5 is typical for dense vegetation
# -----------------------
# minimum_threshold = 2.0
# maximum_threshold = 2.5
# threshold_plot(evi, minimum_threshold, maximum_threshold, width = 8)

# Fractional Cover (pv,npv,bs) = 0 to 100
# Bare Soil (bs) >40 = urbanization / bare soil
# ----------------------
# minimum_threshold = 40.0
# maximum_threshold = 100.0
# threshold_plot(frac_classes.bs, minimum_threshold, maximum_threshold, width = 8)

# Fractional Cover (pv,npv,bs) = 0 to 100
# Vegetation (pv) >80 = dense green vegetation
# ----------------------
# minimum_threshold = 80.0
# maximum_threshold = 100.0
# threshold_plot(frac_classes.pv, minimum_threshold, maximum_threshold, width = 8)
In [30]:
def threshold_count(da, min_threshold, max_threshold, mask = None):
    def count_not_nans(arr):
        return np.count_nonzero(~np.isnan(arr))
    
    in_threshold = np.logical_and( da.values > min_threshold, da.values < max_threshold)
    
    total_non_cloudy = count_not_nans(da.values) if mask is None else np.sum(mask.values)
    
    return dict(total = np.size(da.values),
                total_non_cloudy = total_non_cloudy,
                inside = np.nansum(in_threshold),
                outside = total_non_cloudy - np.nansum(in_threshold)
               )    
    
def threshold_percentage(da, min_threshold, max_threshold, mask = None):
    counts = threshold_count(da, min_threshold, max_threshold, mask = mask)
    return dict(percent_inside_threshold = (counts["inside"]   / counts["total"]) * 100.0,
                percent_outside_threshold = (counts["outside"] / counts["total"]) * 100.0,
                percent_clouds = ( 100.0-counts["total_non_cloudy"] / counts["total"] * 100.0))
In [31]:
# Select a threshold statistical function that matches your land spectral variable
# COUNT = number of pixels in each category
# PERCENTAGE = percent of pixels in each category
# ---------------------------------

# NDBI Threshold
threshold_count(ndbi,minimum_threshold,maximum_threshold, cloud_mask_composite)
# threshold_percentage(ndbi,minimum_threshold,maximum_threshold)

# NDVI Threshold
# threshold_count(ndvi,minimum_threshold,maximum_threshold)
# threshold_percentage(ndvi,minimum_threshold,maximum_threshold)

# EVI Threshold
# threshold_count(evi,minimum_threshold,maximum_threshold)
# threshold_percentage(evi,minimum_threshold,maximum_threshold)

# Fractional Cover - Bare Soil
# threshold_count(frac_classes.bs, minimum_threshold, maximum_threshold)
# threshold_percentage(frac_classes.bs, minimum_threshold, maximum_threshold)

# Fractional Cover - Photosynthetic Vegetation
# threshold_count(frac_classes.pv, minimum_threshold, maximum_threshold)
# threshold_percentage(frac_classes.pv, minimum_threshold, maximum_threshold)
/env/lib/python3.10/site-packages/rasterio/warp.py:344: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned.
  _reproject(
Out[31]:
{'total': 110124, 'total_non_cloudy': 110124, 'inside': 306, 'outside': 109818}
In [32]:
threshold_percentage(ndbi,minimum_threshold,maximum_threshold)
Out[32]:
{'percent_inside_threshold': 0.2778685845047401,
 'percent_outside_threshold': 99.72213141549527,
 'percent_clouds': 0.0}

GeoTIFF Output Products¶

In [33]:
from ceos_utils.data_cube_utilities.dc_utilities import write_geotiff_from_xr

# Remove the comment to create a GeoTIFF output product
# Change the name of the output file, or it will be overwritten for each run 
# Change the desired bands at the end of the function

# Fractional Coverage
# write_geotiff_from_xr("geotiffs/frac_classes.tif", frac_classes, bands=['bs'])

# NDVI
# write_geotiff_from_xr("geotiffs/ndvi_land.tif", ndvi)

# EVI
# write_geotiff_from_xr("geotiffs/evi.tif", evi)

# WOFS
# write_geotiff_from_xr("geotiffs/wofs.tif", water_classification.wofs)
In [34]:
# !ls -lah geotiffs/*.tif